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We report the nonperturbative behavior of the twisted Polyakov loop (TPL) coupling 
constant for the SU(3) gauge theory, which is one of the nonperturbative renormalized 
coupling constants defined in finite volume. We reveal the vacuum structures and the 
phase structure for the lattice gauge theory with the twisted boundary condition, and 
derive the nonperturbative beta function in this renormalization scheme and carry out 
the numerical simulations in the cases of the quenched QCD and Nf ~ 12 SU(3) theories. 

At first, we study the quenched QCD theory by using the plaquette gauge action. 
The TPL coupling constant shows a fake fixed point in the confinement phase even in 
the quenched QCD. We discuss this property and show the nonperturbative running 
coupling constant in the deconfinement phase, where the magnitude of the Polyakov 
loop shows the nonzero values. 

We also investigate the system coupled with fundamental fermions. We use the naive 
staggered fermion in our simulation, and then the minimum number of flavor is 12 in 
this lattice setup because of the twisted boundary condition. The Nf = 12 SU(3) gauge 
theory is expected that the running coupling constant shows a conformal behavior from 
the perturbative analysis. We show the vacuum structure and phase structure of this 
lattice setup with analytical and numerical methods, and then show the detailed analysis 
of the running coupling constant of this theory. Finally, we find the infrared fixed point 
(IRFP) and discuss the robustness of the nontrivial IRFP of many flavor system under 
the change of the analysis method. 

A part of preliminary results was reported in the proceedings [1, 2] and the letter 
paper [3]. In this paper we include a review of these results and show the final conclusion 
for the existence of IRFP of SU(3) Nf = 12 massless theory using the updated data. 

Subject Index B32, B38, B44 



1. Introduction 

Lattice gauge theory is one of the most successful regularization tools to understand QCD. It 
can be applied even to the strong coupling region, and we can investigate the nonperturbative 
properties of QCD. Based on the success of the lattice QCD, there are several application 
for the lattice gauge theory with the different gauge group, the different number of flavor 
and the fermion representation. Among recent lattice studies, the search for the conformal 
or nearly conformal field theory in the infrared (IR) regime have been performed, motivated 
by both theoretical and phenomenological interests [4] - [49]. 

A simple way to find the conformal fixed point is to study the beta function of the theory. 
While the renormalized coupling constant itself is scheme dependent quantity, the existence 
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of the fixed point does not depend on the renormahzation scheme. This is understood as 
follows. Let us consider a relationship between two renormahzed coupling constant defined 
in different renormahzation schemes, gi = F[g2). The beta functions of these two schemes 
are related 

= |£/3(52). (1.1) 

A zero point of the beta function is thus scheme independent except for singular point of 
the transformations. 

In the recent studies to search for the infrared fixed point (IRFP), in particular, there are 
many independent studies in the case of SU(3) gauge theory coupled with Nf = 12 fundamen- 
tal fermions. The existence of the IRFP in Nj = 12 theory was predicted by the perturbative 
2-loop beta function [50] and the further higher order analysis in the MS scheme [51]. The 
phase structure at IRFP was also studied in the paper [52]. Based on these analytical studies, 
there are several works to investigate the nonperturbative running coupling constant and 
the phase structures, although the results have been controversial hitherto. In Ref [8], the 
running coupling constant was computed in the Schrodinger functional (SF) scheme [53-55], 
and within error exhibited scale independent behavior in the IR at coupling g'^'p ~ 5. And 
the MCRG method[22, 23], the phase structure in the finite temperature system [17, 18] and 
the scaling behavior with mass deformed theory [10, 27] show the evidence of the IRFP. In 
the studies of SU(3) gauge theory with Nj = 7 and 10 [4, 41], they found that theory is in 
the conformal window which also suggest that Nf = 12 theory is conformal. On the other 
hand, the studies of the mass scaling behavior [15] and the spectrum of the Dirac operator 
and the chiral symmetry [12, 21] show evidence that this theory is not conformal at low 
energy. 

One reason of the controversy may be the systematic error. In the Ref. [8], they derived the 
running coupling constant by taking the constant continuum extrapolation. That means the 
discretization effects, which is the renormahzation scheme dependent, is neglected. There 
is no study of the running coupling constant which takes care of the discretization error 
carefully at least in the case of Nj = 12. The other reason may be the effect the bulk phase. 
It also depends on the lattice action. Actually, a new bulk phase in this lattice gauge theory 
and the importance of the avoiding bulk phase are also reported [19-21, 24]. 

According to these studies, it is suggested that there are two important things to verify 
the existence or non-existence of IRFP. The first one is to estimate the discretization effect, 
since the strong coupling region might have the large lattice artifact. The critical behavior 
appears in the continuum theory, and then careful continuum extrapolation is necessary. 
The second one is to study the parameter region, in particular to avoid the lattice artifact 
phase. The parameter region of the artifact phase depends on the lattice action, it should 
be studied for each lattice setup. 

This work reports a study of the phase structure and the running coupling constant for 
SU(3) gauge theory with Nj = and 12. We use the plaquette gauge and the naive staggered 
fermion actions. Firstly, we study the phase structure of these theories within both analytical 
and numerical methods, and then derive numerically the running coupling constant in the 
deconfinement phase . The renormahzation scheme we chose is the twisted Polyakov loop 
(TPL) scheme. The TPL coupling was proposed by de Divitiis et al. [56] for the SU(2) case, 
and we extend it to the SU(3) theory. This renormahzation scheme has no 0{a) discretization 
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error, which is of great advantage when we take the continuum hmit. Another advantage 
of this scheme is the absence of zero mode contributions thanks to the twisted boundary 
condition [57, 58]. This regulates the fermion determinant in the massless hmit, which enables 
simulation with massless fermions. In this work, we take the continuum limit carefully, and 
show the existence of the IRFP in the Nf = 12 theory if we include the systematic uncertainty 
coming from the continuum extrapolation. 

In this paper, we give a definition of SU(3) TPL renormalized coupling, and the tree 
level calculation of the TPL coupling in Sec. 2. We show the first example of the running 
coupling constant in the case of quenched QCD theory in Sec. 3. We show the scaling 
behavior of the scheme and the nonperturbative property of the running coupling constant. 
The TPL coupling constant approaches to a constant if the theory is in the confinement 
phase. We discuss the way how to distinguish such a fake "fixed point" and the true one in 
this renormalization scheme in the Sec. 3.1. In Sec. 4, we discuss the vacuum structure and 
Z]y^ symmetry in the system coupled with fermions by the semi-classical analysis in the case 
of Nf = 12 theory. We also show the numerical results of the phase structure of massive and 
massless Nf = 12 SU(3) theory in Sec. 6. In Sec. 7, we study the running coupling constant 
for massless A'^^ = 12 case. We firstly study the global behavior of the step scaling function 
using the data set given in Appendix A, and then investigate the existence of the IRFP by 
the local fit using the additional data set (Appendix B) in the strong coupling regime. The 
detailed discussion for the stability of the IRFP is given in Appendix C. We also discuss the 
taste breaking effects in this analysis in Appendix F. 

One of the main results of this paper is that we found the IRFP at 

5fpL = 2.69 ± 0.14 (stat.)l° 16 (syst-), (1-2) 
and the critical exponent of the (3 function around the IRFP is 

7g* = 0.57t[if (stat.)to,6(syst.)- (1-3) 

There is an independent paper using the similar idea [59]. The present work is that we add 
the new data showing the strong evidence of the IRFP beyond the systematic uncertainty. 
Data analysis is also refined in various ways as discussed below. The details of difference 
from the paper [59], including the discrepancy of the value of ^xpl more than 2-a, are 
discussed in Appendix D and E. 



2. Twisted Polyakov loop (TPL) scheme 

One of nonperturbative definitions for the renormalized coupling constant can be given by 
a divergence-free ratio (Anp) of nonperturbative amplitudes. If the tree level value of the 
quantity is proportional to the squared bare coupling constant, Atree = kgQ, where k is the 
constant which is calculated by the tree level quantity, then we can define the nonperturba- 
tive renormalized coupling constant from the nonperturbative ratio ^atp by identifying the 
renormalization factor of the amplitude as the quantum correction of the coupling constant: 



2 _ Anp 



9np ^ (2.1) 

Since the lattice simulation gives us the value of Anp, what we have to do is to find a ratio 
of tree level amplitudes Atree which is proportional to the squared bare coupling constant. 



3/41 



Twisted Polyakov loop (TPL) scheme is one of such nonperturbative renormahzed cou- 
phngs defined in finite volume. This scheme is given in Ref. [56] in the case of SU(2) gauge 
theory, choosing the ratio of Polyakov loop expectation values for twisted and untwisted 
directions as the quantity Anp- We extend the definition in Ref. [56] to the SU(3) case. This 
scheme can be defined also in the continuum finite volume, however, in this section, we start 
a brief review of the definition of TPL scheme on the lattice. 

2.1. The definition of TPL scheme in the SU(3) gauge theory 

To define the TPL scheme, we introduce twisted boundary condition for the link variables 
(^7^) in X and y directions and the ordinary periodic boundary condition in z and t directions 
on the lattice: 



U^{x + OL/a) 

for fj. = x,y, z,t and u = x,y. Here, {v 
following properties: 

and 



= n,u^{x)nl (2.2) 

= X, y) are the twist matrices which have the 
= I,Tr[J7^] =0, 



for a given fi and u{y^ fi). The gauge transformation f7^(r) — )• A(r)C/^(r)A^'(r + /i) and 
eq.(2.2) imply 

A{r + i)L/a) = nuAir)nl. (2.4) 

In the system coupled with fermions, we also have to define the twisted boundary conditions 
for fermions. The naive twisted boundary condition for lattice fundamental fermions can be 
written by 

ip{x + uL/a) = ni,'ilj{x), (2.5) 

for u = x,y. However, this results in an inconsistency when changing the order of 
translations, namely, 

ip{x + vL/a + pL/a) = rLpVty'ilj{x), 

/ J]^f)p^(x), (2.6) 

for p,v = x,y. To avoid this difficulty, we introduce a "smell" symmetry [61], which can 
be realized by an integral multiple of color symmetry, Ng = Ug x Nc, where Nc and Ng are 
the degree of color and smell symmetry respectively. We identify the fermion field as a 
Nc X Ng matrix (-0^(3;)), here a and a denote the indices of the color and smell and they 
run a = 1, • • • ,Nc and a = 1, • • • ,Ns. Then we impose the twisted boundary condition for 
fermion fields as 

C(x + i>L/a) = e'^/'nf4{n,)l^ (2.7) 

for V = x,y directions. Here, the smell index can be considered as a "flavor" index, then the 
number of flavors should be a multiple of Ng{= Nc = 3). We use staggered fermion in our 
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simulation. This contains four tastes for each flavor. This restrict the number of flavors to 
multiples of 12 in this SU(3) gauge theory with twisted boundary condition. 

Finally, we define the renormalized coupling constant of the TPL scheme. Because of this 
twisted boundary condition, the definition of Polyakov loops in the twisted directions are 
modified as. 



P,(y,z,t) = Tr [l[a,{x=j,y,z,m,e'^^y/^'^ , 



(2. 



in order to satisfy gauge invariance and translational invariance. The renormalized coupling 
in the TPL scheme is defined by taking a ratio of Polyakov loop correlators in the twisted 
{Px) and untwisted (P^) directions: 



2 1 {Zy,,PAy,z,L/2a)Px{o,o,o)^) 

3^PL hatt {Ex,y Pzix, y, L/2a)P,(0, 0, 0)t) ' 



(2.9) 



At tree level, this ratio of Polyakov loops is proportional to the bare coupling. The 
proportionality factor on the lattice (kiatt) is obtained by analytically calculating the one- 
gluon-exchange diagram. To perform this analytic calculation, we choose the explicit form 
of the twist matrices [62] , 
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The Feynman rule for the SU(A''c) gauge theory on the lattice with the twisted boundary 
condition is given in Appendix B in the paper [56] . The value of kiatt is given as 



h 



1 1 



att 



E 



exp(iA;P^ 



<72iV,L2^ V sin2(fc^/2) 



(2.11) 



here L = L/a^ r = {x,y,z,L/2a) and denotes the momentum in each direction. In the 
twisted direction, it gives the sum of the physical and the unphysical twisted momentum: 



k 



x,y 



x,y ' x,yi 



2-Kn 



L 

27rn' 



P% ^ 7r{2mly + 1) 



ph 
z,t 



L 



3L 



(2.12) 



where n^^ = 0, • • • L/2 — 1 and rrix^y = 0, 1, • • • ,Nc — I with {rrix^my) ^ (0, 0). The momen- 
tum k^ can be identified as the color degree of freedom (A'^^ — 1) in the color basis (see: the 
Appendix B in the paper [56]). 
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L/a 



hatt 



4 



0.03213022128143844 
0.03196454161502177 
0.03191145402091543 
0.03188777626443608 
0.03187515361346823 
0.03186277699696222 
0.03185710526062057 



6 



8 



10 



12 



16 



20 



Table 1 The value of hatt for L/a = 4, 6, 8, 10, 12, 16, 20 for SU(3) gauge theory. 

In the continuum limit the proportional factor (fc) for SU(3) gauge theory is calculated 
analytically as: 



The values of kiatt in Table 1 can be fitted by a linear function of O(a^) instead of 0(a), as 
expected. 

3. The TPL coupling for the quenched QCD 

In this section, we study the TPL coupling constant in the quenched QCD. In the first 
subsection, we discuss the property of TPL coupling in the quenched theory and investigate 
the parameter region where the study of the TPL coupling makes sense. We show the running 
coupling constant for the quenched QCD in Sec. 3.2. 

3.1. Phase structure and TPL coupling constant 

The TPL coupling constant is defined by taking the ratio of the Polyakov loop correlators 
between the twisted and the untwisted directions. If the theory is in the confinement phase, 
for sufficiently the correlation length of the Polyakov loop is shorter than the volume and 
the gluon does not feel the boundary effect. We can expect that the ratio of the Polyakov 
loop correlators becomes 1, and we can not extract of the running behavior in such region. 
The quenched QCD theory shows the confinement /deconfinement phase transition in the 
finite volumes, and we can use the TPL running coupling only in the deconfinement phase, 
where the magnitude of the Polyakov loop shows nonzero values. 

To see this property, we observe /3 dependence {(3 = Q/g^ where go is the bare coupling 
constant) of the TPL coupling constant at fixed lattice sizes. Apart from discretization 
errors the scale increases as the /3 decreases at a fixed lattice size. In this test, we use smaller 
lattice sizes, L/a = 2 - 6, with relatively low /3 values. The configurations are generated by 
the Hybrid Monte Carlo algorithm and the action is the Wilson plaquette. We measure the 
Polyakov loop and its correlator for every Monte Carlo trajectory, and each data statistics 
has the same: 20, 000 trajectories. The size of error bar shows the statistical fiuctuation at 
the same statistics. 



k 





= 0.03184- •• . 



(2.13) 



6/41 



The measured TPL coupling constant and the absolute value of the Polyakov loop in t- 
direction are presented in Fig. 1^. The top panels denote the absolute values of the Polyakov 
loop and the bottom ones denote the corresponding TPL coupling scaled by the coefficient 
kiatt for each lattice size. We found that the absolute value of the Polyakov loop approaches 
zero in the low energy region. The confinement /deconfinement phase transition occurs and 
the transition point of /3 depends on the lattice size. From the bottom panels, we can see 
the ratio of Polyakov loop {kiattO^PL) becomes unity below the transition point. 



L/a=2 L/a=4 L/a=6 




2 4 6 85678956789 

P P P 

Fig. 1 The ratio of Polyakov loop and the absolute value of Polyakov loop in t direction 
for L/a = 2, 4 and 6 . 

There is a question whether we can use this TPL scheme for the conformal fixed point 
search in IR region. One discrimination method is to check the the value of renormalized 
coupling, although the value of the renormalized coupling at the fixed point is scheme depen- 
dent. Assume that a theory has IRFP. The fake fixed point appears at A'tpl ~ 1/^ ~ 32. If 
there is an IRFP at ^xpl ~ 32, then we can discuss the fixed point as a physical fixed 

point. On the other hand, if the true fixed point value of the renormalized coupling in the 
TPL scheme is larger than then we will not reach that point. The other important 
check is to see the phase structure of the theory at the same time. At the true conformal 
fixed point, the theory must be in the deconfinement phase. There is a possibility of the 



^ We drop the data of the ratio of Polyakov loop for /3 = 1.0, L/a ~ 2, /3 = 5.0, L/a ^ A and /3 < 
5.5, L/a = 6, since these error bars become huge. Of course, they are consistent with 1. 
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6 7 8 9 10 

P 

Fig. 2 TPL renormalized coupling in the each f3 and L/a in quenched QCD. Dashed lines 
express fit lines of fixed L/a as a function of /3. 

existence of the bulk phase in the low /3 region, in which the Polyakov loop shows the con- 
finement and/or chiral symmetry broken behavior [19, 20, 24], if the lattice theory includes 
the dynamical fermions. We discuss this point in the case of Nf = 12 in Sec. 6. 

3.2. Running coupling constant for quenched QCD 

Now, we would like to show the running coupling constant for quenched QCD. The main 
result was already reported in [1] Here we would entend discussions. The gauge configura- 
tions in this subsection are generated by the pseudo-heatbath algorithm and overrelaxation 
algorithm mixed in the ratio 1:5. One combination of the pseudo-heatbath and 5 overrelax- 
ation steps called a "sweep" in the following. In order to generate the configurations with 
the twisted boundary condition we use the trick [58] proposed by Liischer and Weisz. To 
reduce large statistical fluctuation of the TPL coupling constant, as reported in Ref. [63], 
we measure Polyakov loops at every Monte Carlo sweep and perform a jackknife analysis 
with large bin size, typically of O(IO^). The simulations are carried out with several lattice 
sizes {L/a = 4, 6, 8, 10, 12, 14, 16) at more than twenty /? values in the range 6.2 < /3 < 16. 
We generate 200,000-400,000 sweeps for each parameter set {13, L/a). 

To investigate the evolution of the renormalized running coupling, we use the step scaling 
method [53]. Firstly we chose a value of the renormalized coupling u=g^p^{l3,a/L) at the 
energy scale /x = 1/L. For each L/a in the set of reference lattice size, we flnd the value of 
/? which produces a given value of the renormalized coupling, u. Then, we measure the step 
scaling function on the lattice 

a/L; s)=glp^{f3, a/sL)|g2^^(^^„/i)=„, (3.1) 

by using tuned value of f3 for each lattice size sL/a. Here, s is the step-scaling parameter. The 
step-scaling function in the continuum limit, a{s,u), is obtained by taking the continuum 
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extrapolation of T,{u, a/L; s): 



a{s,u) = limS(n, a/L; s)|g2^j/3,a/i,)=„. 



(3.2) 



This step scaling function {a{s,u)) corresponds to the renormalized coupling at the scale 
fi = 1/sL. 

Figure 2 shows the /3 dependence of the coupling constant in TPL scheme at various lattice 
sizes. The results are fitted at each fixed lattice size to the interpolating function which is 
similar to the one used in Ref . [8] , 



tip' 



B) 



(3.3) 



where Ai are the fit parameters, and 4<i?<5, n = 3, 4 are employed. As reference lattice 
sizes of the step scaling, we use L/a = 4, 6, 8, 10. The step scaling parameter is s = 1.5, and 
we estimate the coupling constant for L/a = 9, 15 from interpolations at the fixed (3 using 
the above fit results of all the lattice sizes. 

We take the continuum limit using a linear function of (a/L)^, because the TPL scheme 
involves no 0{a/L) error. We found that the coupling constant in the TPL scheme exhibits 
scaling behavior even at the smaller lattice sizes, as shown in Fig. 3. 
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Fig. 3 The continuum extrapolations of ^tpl scale L/Lq = s". The fit function is 

a linear function of (a/L)'^. The data from left to right correspond to L/a = 15, 12, 9, 6. The 
statistical error bars are of the same size of the symbols. 



The running of the TPL coupling constant in quenched QCD with 24 steps is shown in 
Fig. 4 together with one- and two-loop perturbative results. The horizontal axis corresponds 
to the energy scale. All the results are normalized at L = Lq with g'^{Lo) = 0.65. The non- 
perturbative running coupling constant is consistent with one- and two-loop perturbative 
results in the high energy region {Lq/L > 0.1). Comparison with the other nonperturbative 
analyses in SF scheme and Wilson loop scheme are also interesting. The both SF (Fig. 1 in 
the paper [64]) and Wilson loop (Fig. 8 in the paper [65] running coupling constants show 
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that similar behavior; faster than one-loop in the nonperturbative region. On the other hand 
the TPL running is slightly slower than one- loop. 




I I I I I I 

0.0001 0.001 0.01 0.1 1 10 

L„/L ill/ A) 



Fig. 4 The running coupling constants in TPL scheme, and in perturbative one-loop and 
two- loop calculations. 



From this quenched test, we conclude that we can control both the the statistical and 
systematic errors of the TPL coupling constant. Furthermore we find that the TPL coupling 
constant in quenched QCD has a robust scaling behavior even in a small lattice size, which 
was also observed in the previous quenched SU(2) calculations [56, 63]. 

4. Vacuum structure of the Nf = 12 theory with twisted boundary condition 

Now, we would like to consider the SU(3) theory coupled to Nf = 12 fundamental fermions. 
In this section, we discuss the center symmetry of SU(3) gauge group to define the true vac- 
uum in this theory. The generators of the center symmetry of SU(A'^c) pure gauge theory are 
z = exp(27ril/Nc), where / = 0, 1, • • • ,Nc — 1- This symmetry is broken by adding fermions 
to the theory, leading to the existence of a unique vacuum in an SU(3) gauge theory involv- 
ing massless fermions in the deconfinement phase. We discuss the vacuum structure which 
is an important to study of the gauge theories in finite volume. 

Let us focus on the center symmetry of this theory. Although the Wilson gauge action is 
invariant under the following transformation for the link variable for each direction, 

U^{t,x) ^ zUf,{t,x), (4.1) 

the fermion is not invariant. 

At the perturbative one-loop level, the semi-classical free energy for the gauge configuration 
{[/} is given as 

^(tree and one-loop) ^ S^lf/) + S'°''''"^°°P[f/] - TV/ In det [/;([/)]. (4.2) 

With the twisted boundary condition, the flat potential due to the toron contribution is lifted 
because of the nonzero momenta in twisted directions and the free energy has 3"^ = 81 -fold 
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-lndet(D) - So 


(0,0) 





(0, 1 or 2) or (1 or 2,0) 


-57.19 


(1 or 2, 1 or 2) 


-89.56 



Table 2 The 6z,t dependence of semi-classical free energy for 6^ lattice. We take a reference 
potential energy Sq = —3511.68 which is the one in {0z,6t) = (0,0). 

degenerate classical minima at Ufj, = exp(27rz^^/3)I, where 6^ = 0,1,2 for each direction. 
The Wilson gauge action (Sg) and the one-loop contribution from gauge part {Sg^^ ^'~"~'^[U]) 
respects the Z3 symmetry, so that we do not consider them in what follows. 

Let us consider the fermion determinant. In the momentum space, there are the physical 
and unphysical momenta {k^ = + k^) in the twisted directions, that also appear in the 
gauge field momenta (eq. (2.12)). In the case of the fermion field, the color and smell degree 
of freedom of in eq.(2.7) can be transferred into the unphysical momentum degrees of 
freedom: their number is A'^c x A'^s ^. Here, we replace the momentum as 

k^^k1^ = k^ + 2Tri9^/3L. (4.3) 

The Z3 transformation in eq. (4.1) can be defined on each lattice site independently, so that 
we can take a typical gauge in which = exp(27ri0^/3L) for whole lattice volume. Then the 
fermion action in the vacuum = exp(27ri0^/3-L) is obtained by the above replacement. 
The fermion determinant in finite volume is thus 

-lndet[L>] = -8^1n [ ^sin2(^^) j . (4.4) 

In the Table 2, we show the results of the fermion determinant (eq.(4.4)) for L/a = 6. We 
find that the vacuum free energy is independent of 6x^y and there are three types of vacuum 
classified with O^^t- The first one is a "trivial vacuum", in which vacuum {Oz,Ot) = (0,0). 
This vacuum has 9— fold degeneracies. The value of the free energy is highest among three 
types of vacuum, so that it will decay to the true vacuum. The second one is a "half-trivial 
vacuum" , in which one of O^^t is 1 or 2 and the other one is 0^ = 0. This vacuum has 36— fold 
degeneracies. The free energy is higher than the one of the third vacuum, so that this vacuum 
is also unstable. The third one, in which the free energy is lowest, is a "non-trivial vacuum". 
Both 6z and 9t take 1 or 2, and there are also 36— fold degeneracies. The free energy has 
minima at this vacuum and in this non-trivial vacuum all classical link variables for z and t 
directions has a non-trivial phase Uz,t oc exp(±27ri/3L). This means that the Polyakov loop 
for z direction also has a non-trivial phase exp(±27ri/3). 

This classification is hold for generic lattice size, and it turns out that the potential def- 
erences between the true non-trivial vacuum and the other vacua become small when the 
lattice size become large, although the potential barrier must become high. If we change 



^ In the case of Ng = rig x iVc, there are Ug flavor fermions whose momentum in the twisted 
directions have x Nf. unphysical modes 
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the fermion boundary condition in z and t direction, then the semi-classical free energy has 
minima at the other vacuum. 

5. Simulation setup for the N ^ = 12 theory 

Our numerical simulation is performed in the following setup. The gauge configurations are 
generated by the Hybrid Monte Carlo algorithm, and we use the Wilson gauge and the 
naive staggered fermion actions. In Sec. 6, the simulations are carried out with lattice sizes 
Lja = 4,8 and 12 at several low (3 and a broad range of ma to study the phase structure 
in this system. We generate 1, 000 - 2, 000 trajectories for each parameter set in the case of 
L/a = 4 and 8, and also generate 500 ~ 1, 000 trajectories for each in the case oi L/a = 12. 

The measurement of the coupling constant in Sec. 7 are carried out with lattice sizes 
L/a = 6,8, 10, 12, 16 and 20 ^ at around thirty /3 values in the range 4.0 < /3 < 100. To 
reduce statistical fluctuations, we measure the Polyakov loops at every trajectory and bin 
the data by taking the autocorrelation into account. Using the jackknife method, typical 
statistical errors of correlator are 2 - 3%. We also estimate the statistical error within the 
bootstrap method in the whole analysis in Sec. 7, whose results are consistent with the 
jackknife analysis. 

6. Simulation results: Phase structure of Nf = 12 SU(3) theory with the 
twisted boundary condition 

In this section, we investigate the phase structure of the Nf = 12 fermion theory on the 
lattice. Although the main purpose of this paper is to derive the running coupling behavior 
in the massless Nf = 12 theory, in order to fully understand the phase structure we need to 
understand the phase structure of the whole region of the theory space including the mass 
parameter. Actually, there are several studies which reported the existence and absence of 
the bulk phases in the case of Nf = 12 staggered fermion system [19, 24]. The paper [24] 
suggested a shift symmetry broken phase in the case of the improved staggered action, and 
the paper [19] reported that there is no such phase in the case of the naive staggered action 
and the shift broken phase comes from the improvement of the staggered fermion action. 
Furthermore it is suggested that there is the spontaneous chiral symmetry broken phase in 
the strong coupling limit for Nf < 52 [20, 66]. In our simulation, we use the naive staggered 
and introduce the twisted boundary. It is important to show the phase structure in our 
lattice setup independently. We identify the parameter range to study the TPL coupling 
constant. 

At first the plaquette values on the /3 — ma plane be reported in Sec. 6.1. In Sec. 6.2 we 
discuss the complex phase of the Polyakov loops and the phase structure of the theory. In 
Sec. 6.3, the tunneling behavior between the degenerate vacua and whose effect to the TPL 
coupling constant are discussed. 



We generated the data with L/a = 4 as same as the quenched QCD case, however, we found there 
is a large discretization effects in the case of Nf = 12 [2], and the systematic uncertainty could not 
be controlled. Then we dropped the data from this analysis. 
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6.1. Plaquette values on the (3 — ma plane 

Let us investigate the plaquette values of the f3-ma plane. The left panel of figure 5 shows 
the plaquette values on (L/a)'^ = 4^ lattice in the fi-ma plane in the range of < ma < 0.2. 
Most of the configurations are thermalized from massive to massless direction except for the 
small mass region in the /5 = 3.8 as we explain later. 
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Fig. 5 Plaquette values for several /3 and ma on (L/a)^ = 4^ lattice. The data at ma = 1 
on the right panel denotes the quenched QCD. The small (red) arrows near the massless at 
f3 = 3.8 on the left panel shows the detailed history of the thermalization. The other data 
are thermalized from massive to massless direction. 



First of all, we found there are gaps of the plaquette in the lines of fixed ma. The gap 
appears between /? = 4.8 and /? = 4.6 for ma = 0.2 and in the smaller mass region it happens 
at the lower /?. Near the massless limit the gap exists around 3.6 < /3 < 4.0. 

The small (red) arrows near the massless at /3 = 3.8 on the left panel in Fig. 5 shows 
the detailed histories of the thermalization. We find that there are two different values in 
the < am < 0.0125 region. The configurations giving larger values of the plaquette at the 
same ma are thermalized from the massless configuration in /3 = 4.0; on the other hand 
those giving smaller values are obtained from the massive configuration at fixed f3. It clearly 
shows the first order phase transition around this region. At /3 = 4.0 and /3 = 3.6, there is 
no dependence on the thermalization process. 

We also study the larger mass region. The right panel in Fig. 5 shows the same plot with 
the left for a broader region of ma. In the quenched limit, we know that there is the first 
order phase transition. In the figure, we plot the data for the quenched lattice at ma = 1.0. 
The gap seems milder in the larger mass region. 



13/41 



L/a=8 



L/a=12 




0.2 0.4 0.6 0.^ 

ma 



beta=3.6 
beta=3.8 
beta=4.0 
beta=4.2 
beta=4.4 
beta=4.6 
beta=4.8 



beta=5.0 

beta=5.25 

beta=5.5 

beta=5.75 

beta=6.0 

beta=6.5 



Fig. 6 Plaquette values for several /3 and ma on (L/a)'^ = 8'^, lattice. The data at 
ma = 1 denotes quenched QCD. The statistical error is the same size of the symbol. 



We also investigate this first order phase transition near the massless region by changing 
the lattice volume in Fig. 6. There are slight differences between L/a = 4 and L/a = 8 for 
the critical value of (3, for example at ma = 0.175 the data for [3 = 4.6 on (L/a)"^ = 4^ is 
in the upper phase of the gap, but it moves to the lower phase in the case of (L/a)^ = 8^. 
The results on (L/a)^ = 8^ and 12^ show the similar behavior at least the present interval of 
P and ma (A/3 = 0.2, Ama = 0.01 - 0.025), and there is no clear volume dependence. The 
massless simulation at /3 = 3.8 needs extremely finer molecular-dynamics time step size than 
At = 0.002 (r = 1 is 1 trajectory), and then practically we could not generate the data. The 
position of /3 where the simulation becomes quite costly is the same in the both L/a = 8 
and L/a = 12. This supports the understanding that there is a bulk phase transition. 



6.2. Polyakov loop 

Next, let us investigate the Polyakov loop. Now, there are dynamical fermions and the 
center symmetry is broken explicitly, and thus there is no clear order parameter for the 
deconfinement phase transition. However, here we use the word "deconfinement" phase as 
where a magnitude of Polyakov loop is clearly nonzero on the lattice. At first, we show the 
complex phase of Polyakov loop. 

In the case of the quenched QCD, there is Z-^ symmetry, and there is tunneling between 
them on finite lattices. In the case of Nf = 12 massless theory with the twisted boundary 
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condition, the symmetry is broken and the true vacua is the one that the Polyakov loops 
in the untwisted directions have the nontrivial phase as explained in Sec. 4. 

We present the scatter plots of the typical Polyakov loops for the massless theory in the 
twisted and untwisted directions in the left and right panels in Fig. 7 respectively. The 




Fig. 7 Scatter plots of the Polyakov loop in the twisted {fi = x) and the untwisted (fi = t) 
directions in the complex plane. The lattice size is 6^ and circle(blue), square(red) and 
cross(green) symbols denote f3 = 20.0, 8.0, 4.0 respectively. 



Polyakov loop in the twisted direction clearly shows the independence of the complex phase. 
That is consistent with the tree level analysis, in which (Px) = because of the traceless 
twist matrix in the eq.(2.8). This does not depend on the value of /3 and it is not related to 
whether the system is in the confinement or in the deconfinement phase. On the other hand, 
the Polyakov loop expectation value in the untwisted direction has clearly the nontrivial 
phase exp(±27ri/3) in the high f3 region. In /3 = 4.0 the tunneling occurs between the two 
phases, but apart from the tunneling the values of the phase are close to exp(ib27rz/3). Thus, 
we confirm that the Polyakov loops in both the directions are consistent with the results of 
from the semi-classical analysis in Sec. 4 even in the strong coupling region. The effect of 
the tunneling on the TPL coupling will be discussed in the next subsection. 

Prom these scatter plots, the real part of the Polyakov loop in t-direction is a negative 
nonzero value in the deconfinement phase because the nontrivial phase vacua are chosen 
as shown in Sec. 3. The left two panels in Fig. 8 shows the real part of the Polyakov loop 
in t-direction at several fixed ma in the case of L/a = 4. The Polyakov loop in the high 
/3 region clearly deviates with zero for all ma, which implies that the system stays in the 
deconfinement phase. On the other hand, at the low /? region that is consistent with zero 
and the theory shows the confinement behavior. In the case of massless L/a = 4 we find a 
gap at the /3 = 3.8, which is at the same position for the plaquette study. For /3 smaller than 
the gap position, the real part of the Polyakov loop is not consistent with zero, but it goes 
to zero continuously. In the finite mass region, there is a weak jump, and the gap become 
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Fig. 8 Left panel: real part of the Polyakov loop in t direction in the case of {L/a^ = 4^- 
Clearly there is a gap for each ma. In the case of the massless, the measured data at /3 = 3.8 
are shifted; the one of them shows at /3 = 3.77 which corresponds to the configurations 
at /3 = 3.8, ma = whose plaquette value is the larger one in Fig. 5. The other shows at 
(3 = 3.83 which corresponds to the configurations at the same /3 and ma whose plaquette 
value is the smaller one in Fig. 5. Right panel: the same plot in the case of (L/a)^ = 8^. 

larger in the smaller mass region. The value of the critical /3 in which the data shows the 
jump is the same with the plaquette study. 
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Fig. 9 Left panel: real part of the Polyakov loop in t direction in the case of /3 = 5.0, ma = 
0.4, {L/a)^ = 4^. Right panel: the same plot in the case of the same /3 and ma on (L/a)^ = 8^. 



Let us study the lattice size dependence. We show the real part of the Polyakov loop in 
the case of L/a = 8 in the right panels in Fig. 8. There is no clear jump in the case of 
L/a = 8, but the real part of the Polyakov loop approach to the zero in the low /? region. 
Note that at large mass region, for example /3 = 6.0, ma = 0.4, we can find that the theory 



16/41 



is in the deconfinement phase from the scatter plot of the Polyakov loop on the complex 
plane, although the real part of the Polyakov loop is consistent with zero. At the parameter 
the fermion mass is too heavy, and then there is the tunneling between Z3 vacua as in the 
quenched case. In this analysis, the number of trajectories is the same for each parameter 
set at the fixed mass and then we can naively find such a fake Re Pj = data from the size 
of error bars and can confirm it by the scatter plot. 

In the top panels of Figs. 8, the theory with ma = 0.4, L/a = 8 is in the confinement 
phase for /3 < 5.0, while in the case of L/a = A the theory is in the deconfinement phase 
for 4.8 < p. Actually, figure 9 is the scatter plots showing a clear difference. On the other 
hand, ma = 0.025 shows the transition between /3 = 4.2 and /3 = 4.0, it is the same position 
in the case oi L/a = 4. At least the current interval of /3 and ma, we cannot find the volume 
dependence in the small mass region. The data for ma = 0.4 and 0.2 show the volume 
dependence, and then it seems to be crossover for both the bulk and the finite volume phase 
transitions. In the case of the (nearly) massless fermion, the transition becomes bulk, and 
the theory with the massless fermion lies in the deconfinement phase for (3 > 4.0. The value 
of the critical (3 is the same with the one of the gap of the plaquette. In Figs. 5 and 6, the 
phase for (3 larger than the gap position can be identified as the deconfinement phase and 
that for (3 smaller than the gap position must be the confinement phase. Furthermore, in the 
case oi L/a = 12, we cannot find clear nonzero value of the Polyakov loop in the whole region 
since the current preliminary statistics is small and the Polyakov loop in the low f3 region is 
noisy. The gap of the plaquette and the confinement /deconfinement phase transition seems 
to occur simultaneously, and the there is no difference between L/a = 8 and L/a = 12 of the 
plaquette behavior. At least we found that the position of /3 where the simulation becomes 
quite heavy is the same in the case of L/a = 8. 

PI 

deconfine « 




confine 



ppig quenched 

Fig. 10 The phase structure of Nf = 12 SU(3) theory with naive staggered fermion. The 
dashed line denotes bulk phase transition and the solid line denotes the finite volume phase 
transition. 

The summary of the phase structure and the available region of the TPL coupling for 
the quenched and the massless Nf = 12 QCD is the following. Figure 10 is a sketch of the 
phase structure for the naive staggered Nf = 12 SU(3) theory. In the case of the quenched 
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QCD, the correlation length becomes shorter in the lower (3 region, and there is the finite 
volume phase transition where the theory goes to the confinement phase. In the case of the 
massless Nf = 12 SU(3) theory, there is the similar behavior while the transition seems the 
bulk one at /3 < 4.0. In both cases, the boundary condition is negligible in the confinement 
phase where the correlation length of the Polyakov loop is smaller than the lattice extent, 
and then the TPL coupling becomes the constant. In the study on the running coupling 
constant in TPL scheme, we should focus on only the deconfinement phase on the lattice. 

Finally, we show that our massless configurations to study the running coupling constant 
live in the deconfinement phase (although it might be trivial since the transition seems to be 
the bulk and we concentrate on the parameter region within (3 > 4.0). Figure 11 shows the 
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Fig. 11 Real part of the Polyakov loop in the untwisted (^u = t) direction as a function of 
the lattice size at /3 = 6.0. The statistical error bars are the same size of the symbols. 



real part of the Polyakov loop for the low (3 region for each lattice size. All data points clearly 
show that Re Pj take the non-zero values. This observation indicates that all configurations 
used in our analysis in Sec. 7 are in the deconfinement phase. 

6.3. Tunneling behavior between true vacua 

During the Monte Carlo simulation, the tunneling can occur between the two degenerate 
vacua in each untwisted direction independently. The tunneling behavior is a lattice artifact, 
since the potential barrier between the vacua becomes finite due to the finite lattice spacing. 
In this subsection, we would like to consider how the TPL coupling is disturbed by the 
tunneling behavior. 

The tunneling is expected to occur more frequently in the strong coupling region. We 
observe some tunnelings in low /3 region, although the number of the tunneling is quite 
small, and decreases as /3 increases. For L/a = 6 we observe the tunneling 7 times per 6, 000 
trajectories in the stronger coupling region, /3 = 4.0, and almost no tunneling, 3 times per 
90, 000 trajectories, at (3 = 6.0. 
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A typical example of the tunneling is shown in figure 12. The left panel denotes the history 
of the imaginary part of the Polyakov loop in the untwisted direction at /? = 4.0 on L/a = 6. 
In the figure the sign of the imaginary part is flipped in the 1,400 - 1, 600th trajectories. 

Let us see the effects of the tunneling on the TPL coupling defined by the ratio of the 
Polyakov loop correlators in eq.(2.9). The right panels in Fig. 12 present the histories of the 



ratio of Polyakov loop correlator 




trajectory 

Fig. 12 Left panel: History of the imaginary part of the Polyakov loop in untwisted {fi = t) 
direction at /3 = 4.0 in L/a = 6. Right panels: Histories of the numerator and denominator 
of the TPL coupling eq.(2.9) in the same trajectories of the left. 

numerator and denominator of the TPL coupling obtained from the same configuration as in 
the left panel. During the tunneling, i.e., the 1,400 - 1,600th trajectories, both the results 
seem not to exceed the range of each statistical fluctuation. This means that the tunneling 
does not significantly affect the result of the TPL coupling calculated by the ratio of the 
expectation values for the numerator and denominator. 

Therefore we conclude that effects of the tunneling on the TPL coupling is negligible at 
least in our calculation due to the rare occurrence of the tunneling and the large statistical 
fluctuation of the Polyakov loop correlators. 

7. Simulation results: Step scaling function for A^^ = 12 SU(3) gauge theory 

In this section, the step scaling function in the case of Nf = 12 flavor are derived. We focus on 
the growth rate of the step scaling function as a function of the input renormalized coupling. 
This quantity is a suitable for showing the precise behavior of the running coupling constant, 
and becomes unity when there is a zero in the beta function. At first, we will discuss the 
global behavior of the growth rate from the perturbative to the IR region in Sec. 7.1. The 
nonperturbative running behavior shows the signal of the conformal fixed point in the IR 
region. Then, we focus on the low energy region only and derive again the step scaling 
function by using the data only in the strong coupling region in Sec. 7.2. We discuss the 
stability of the IR fixed point by considering several systematic uncertainties and derive the 
universal quantity for the exponent of the beta function around the IRFP in Sec. 7.3. 
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7.1. The global behavior of the step scaling function 

Before explaining the simulation details, we address the guiding principles of our simulation 
to show the global behavior of the running coupling. We use the step scaling method as 
same as the quenched case. For each L/a we use the interpolation in (3. Since the choice of 
the interpolation function can affect the step scaling function, we should take care of the 
following two points: 

(1) We generate data in a broad range of /3, with intervals that the renormalized coupling 
constant (5^) grows almost constantly in each interval. Thus the interval of /3 is large 
in high /3 region while small in the low /3 region. Each data have a similar accuracy 
(2-3%). 

(2) We employ fit functions for /3 interpolation which reproduce the tree level result gl^ ~ 

on each lattice size in extremely high /3 region. 

These guiding principles ensures the stability of our fit results under the change of fit func- 
tions and the number of data. Point 1 is needed to ensure that the fit functions do not favor 
any special region of the data when we interpolate our data in /3 or extrapolate to in (a/L)^. 
To satisfy this point is very important to search for the IRFP by using the numerical simu- 
lations, since finally the interpolation function of the data decides the step scaling function. 
We discuss the dependence of the data set for the global fit analysis in Appendix D, in the 
case where the data are concentrated in some particular region. The result in Appendix D do 
not have a nice agreement with the perturbative result and in the IR region the position of 
the IRFP strongly depends on the fit range. Point 2 is needed to reduce the effect of statis- 
tical fluctuation in high (3 region. In high energy region, since the coupling runs very slowly, 
we need extreme high statistics to reproduce the perturbative result only by the numerical 
data. The assumption of the point 2 makes the result stable in high energy region. 

In Fig. 13, we show our simulation results for the renormalized coupling in TPL scheme as 
a function of 1//3 for each L/a. The raw data are given in Tabs. Al and A2 in Appendix A. 
The left panel in the Fig. 13 shows a global behavior of the TPL coupling. We can see the 
high energy behavior seems to be almost linear in 1//3 as expected from the perturbation 
theory. The right panel focuses on the low (3 region. In the low (3 region for L/a = Q the 
TPL coupling has a maximum at /3 = 4.3. In the contrast to the Schrodinger functional 
scheme [8], there is no inversion of the lattice size dependence at fixed (3. We consider that 
this difference comes from the lattice artifact which depends on the renormalization scheme. 
To remove the effect, the careful continuum extrapolation is necessary. 

For /3-interpolation, we use the following form of fitting function: 

gW{P,a/L) = 6//3 + Y.ti a(a/L)//3^+\ (7.1) 

where Ci{a/L) are the fitting parameters and N is the degree of the polynomial. Here, = 3 
- 5 are employed. We drop /3 = 4.0, L/a = 6 from the fit to avoid the double solutions when 
we solve the (3 values to reproduce the input renormalized coupling (u). This limits us to 
study the step scaling in the range u < 2.94, where u = 2.94 is the renormalized coupling 
constant ^^PL at /3 = 4.3, L/a = 6, in this lattice set up. 

To investigate the evolution of the renormalized running coupling, we use the step scaling 
method as explained in Sec. 3.2. In this study, we take s = 1.5, and denote cr{u) = a{s=1.5, u) 
in the rest of this paper for simplicity. The set of small lattice is taken to he L/a = 6, 8, 10, 12, 
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Fig. 13 TPL coupling for each /3 and L/a. Left panel: Plots for the global region of /3. 
Right panel: Plots only for the low /3 region. 



therefore, we need values of ^tpl L/a = 9, 12, 15, 18 to take the continuum limit in 
Eq. (3.2). For L/a = 9, 15 and 18, we estimate values of ffxPL ^ given /? by the lin- 
ear interpolation in (a/L) with using the data on the lattices L/a = {8, 10}, {12, 16} and 
{16,20}, respectively. To estimate the systematic error of this interpolations, we also per- 
formed the linear interpolation in (a/L)^, and found that the difference between a/L and 
(a/L)'^ interpolations is negligible. Furthermore, we will also show the s = 2 step scaling in 
which there is no interpolation of L/a, and the difference of the result should give an indirect 
estimation of the systematic error from the interpolation in L/a. 
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Fig. 14 Continuum extrapolation for the case of input coupling u = 0.5, 1.5 and 2.5 from 
the bottom to top respectively. Each solid line denotes the input renormalized coupling u, 
and the dashed line with the corresponding color denotes the linear extrapolation function 
of (a/L)2. 
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In Fig. 14, we show the examples of the continuum extrapolation for obtaining cr{u) in 
the weak, intermediate and strong coupling regions. We determine the central value of ^{u) 
by the linear extrapolation in {a/ L)"^ with four points; L/a = 6, 8, 10, 12 — )■ 9, 12, 15, 18 since 
the leading order of the discretization error comes from O(a^) in this scheme. Note that, in 
the strong coupling region, each lattice data S(u, a/L;s = 1.5) (black data point) is quite 
larger than n, however, at the continuum limit, (t{u) gets close to u. This indicates that it is 
very important to take the continuum limit carefully in this kind of analysis. We explain the 
reason why this continuum extrapolation is chosen as the best in this analysis in Appendix E. 
We also discuss the taste breaking in the continuum limit in Appendix F. 

1.06 



1.04 



1.02 



0.98 



Fig. 15 The growth rate a{u)/u as a function of u with statistical error. Two-loop pertur- 
bative value (black line) is also plotted for comparison. The horizontal (green) line denotes 
unity line, where the beta function is consistent with zero. 

Now, we obtain the step scaling function explained above in a wide range of u. Figure 15 
shows the growth rate of the renormalized coupling [(t{u) /n) as a function of u with statistical 
error which is estimated by jackknife method. We also carried out the bootstrap analysis 
independently, and found that the results are consistent with this jackknife analysis. 

We found two things from this plot. The first one is that the result is consistent with 
perturbation theory in the weak coupling regime. The TPL coupling coupling with this 
lattice set up works promising under this analysis method. The other one is the central 
value of a{u)/u becomes unity around u = 2.7, demonstrating the signal of a fixed point. 
This is the first zero of the beta function from the asymptotically free regime. It suggests 
the existence of an infrared fixed point around the region. Unfortunately, the upper values 
of the error bars do not cross the line a{u)/u = 1. This means that we cannot exclude 
the possibility for the coupling constant to continue growing within the error bar. We will 
investigate this quantity again by focusing only the strong coupling region and adding the 
data. Furthermore, will give an estimation of the systematic error of the fixed point coupling 
in the next subsection. 

7.2. Low energy behavior and stability of the IR fixed point 

In the previous subsection, we found a signal of the IRFP around u = 2.7 from global fit 
of the data. Now we focus on the strong coupling region and will determine the fixed point 
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7^ of the data 
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4.3 < /3 < 7.0 
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Table 3 The fit range of the local fit analysis for each lattice size. The data L = 9, 15 and 
18 are obtained by the L/a interpolation at the fixed j3 value as explained in Sec. 7.1. 



coupling and the related universal quantity. In this subsection, we take a narrow (3 range 
in which /3-dependence of 9xpl approximated by linear or quadratic functions of /?. 

We add more data, which is a part of the data shown in Appendix B, to obtain the precise 
result and discuss the systematic uncertainties of the IRFP. 

Practically, we will carry out the step scaling again with the data only in low (3 region 
u > 2.0. This region roughly corresponds to the range /3 < 7.0. We choose the fit range for 
each lattice size as shown in Table 3. The fitting function is chosen as a simple unconstraint 
polynomial function: 

= E£o' C,{a/L)/I3\ (7.2) 

where is the degree of the polynomial and here we take = 3 for L/a = 6, 8 and N = 2 
for the other lattice size. We derived the step scaling function by using the same procedure 




Fig. 16 The local fit result for the growth rate of the TPL coupling. The solid (blue) error 
bar denotes the statistical error and the dot (black) error includes the systematic error. Two- 
loop perturbative value (black line) is also plotted for comparison. The horizontal (green) 
line denotes unity line, where the beta function is consistent with zero. 



with the pervious subsection. The growth rate of the step scaling function is shown in Fig. 16. 
As a central analysis with solid blue error bar, we take the four point linear extrapolation 
in (a/L)^ with statistical error estimated by the jackknife method. The dot (black) error 
bar includes the systematic error, we will discuss it soon later. This local fit result clearly 
crosses the line a{u)/u = 1, which shows the existence of the IRFP. 
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Fig. 17 The comparison of the results with the global fit and those with the local fit. Two- 
loop perturbative value (black line) is also plotted for comparison. The horizontal (green) 
line denotes unity line, where the beta function is consistent with zero. 

Figure 17 shows the comparison between the local fit result with the global fit result. These 
two central value are consistent with each other within 1-cr, despite the change of the data 
set, the fit range and the fitting function. The results strongly consolidates the existence of 
the stable fixed point around u = 2.7. The error bar for the local fit analysis is smaller since 
owning to the additional precise data. We also report the data sets dependence and the fit 
range dependence independently in Appendix C. 

Now, we would like to estimate the systematic error in our analysis. There are two possible 
dominant sources of the systematic error. One is from the choice of the fit range for the (3- 
interpolation (Eq. (7.2)). As shown in Fig. 17, there is a small difference between the global 
fit and the local fit in Table 3, and we also investigate the further narrow range of the /3. 
Even if we take only the data of /3 < 7.0 for all lattice size, the fitting function does not show 
a large difference. We can conclude the systematic error from the choice of the fit range is 
small in this analysis (see Fig. CI in Appendix C). 

The other dominant systematic error comes from the continuum extrapolation. In Fig. 18, 
we show the comparisons of several types of continuum extrapolation for u = 2.0, 2.686 and 
2.85. As the central value, we take the linear extrapolation in [ajV)^ for Lja = 6, 8, 10, 12 — t- 
9, 12, 15, 18. We estimate the systematic error by taking the difference between the central 
value and the result from linear extrapolation without the data on the coarsest lattice L/a = 
6. Furthermore we compare the central value with the quadratic extrapolation with all the 
data at four L/a. Figure 18 shows the TPL renormalized coupling has a small systematic 
error in the strong coupling region, and all the values in the continuum limit agree within 
1-0" statistical errors. The total error in Fig. 16 is estimated by adding the difference between 
the continuum extrapolations as a systematic error to the statistical error in quadrature. We 
conclude that the existence of the IRFP is stable in this analysis. 

Finally the renormalized coupling at the IRPF is 

u* = 2.69 ± 0.14 (stat.)_o.i6 (syst.). (7.3) 



2-loop 

■ — ■ global fit 
•— • local fit 
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Fig. 18 Continuum extrapolation for the case of input coupling u = 2.0, 2.686 and 2.85. 
Each green line and the blue curve denote the 4 points linear and quadratic extrapolation 
functions of O(a^) respectively. The red line show the linear extrapolation function of 0{o?) 
for 3 data points without the coarsest lattice data. In the case oi u = 2.0, the step scaling 
function is larger than the input value, however, it becomes consistent with u aX, u = 2.686 
and for the larger u it is smaller than the input renormalized coupling constant. 



Here, the jackknife error of the running coupling is used as a statistical error and we estimated 
the systematic error coming from the continuum extrapolation. The corresponding /3 value 
for each lattice size at u* = 2.69 can be calculated from the (3 interpolation function in 
Eq. (7.1): (/3, L/a)=(4.91,6), (5.41,8), (5.65,10), (5.79,12), (5.91,16), (5.94,18), (5.94,20). 
These are the parameter sets which reproduce conformal physics after taking the continuum 
limit. 

In addition, we mention the numerical stability of our results. For this purpose, we perform 
another step-scaling analysis based on s = 2 with L/a = 6, 8, 10 — )• 12, 16, 20. The continuum 
limit is taken by linearly extrapolating three points in (a/L)^. The explicit figure of the 
growth rate of the renormalized coupling is shown in Fig. CI in Appendix C. We find that 
the running behavior in s = 2 step scaling is also consistent with that in the case of s = 1.5, 
and the IRFP is found at u* = 2.49 it 0.19 (stat.). This is consolidating the existence of the 
IRFP in this theory. 

7.3. Critical exponent 

Finally we will derive the critical exponent around the IRFP. In this theory, we have one 
irrelevant parameter, which is the renormalized coupling constant, around the nontrivial 
fixed point. The value of the renormalized coupling at the fixed point depends on the renor- 
malization scheme. If we denote the scheme transformation from one renormalization scheme 
ui to the other one U2 = F{ui), then in the perturbative region, the function F{ui) can be 
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expanded as a polynomial function. The beta function for the renormalized coupling is uni- 
versal up to two-loop order, however, the nonperturbatively the one for U2 scheme is related 
with the other one as follows: 

/5(n2) = (7.4) 

In the vicinity of the IRFP, the beta function can be approximated by 

/?(n) ~ -7g*K -u) + 0{{u* - uf). (7.5) 

Although the value of renormalized coupling at the IRFP is scheme dependent, we can easily 
find the coefficient 7* around the point is the scheme independent quantity using eq. (7.4). 

Now, we compute 7* from the slope of a{u)/u against u, and obtain s~^!i = 0.79 it 
0.11(stat.) in the central analysis in the Fig. 16. This leads to 

% = 0.57l[jJ5(stat.)_o.i6 (syst.), (7.6) 

where the first error is statistical error using the jackknife method and the second one is 
the systematic error from the continuum extrapolation estimated by the comparison to the 
3 point linear continuum extrapolation. The value of 7* is sensitive to the variation of the 
slope, which causes rather large statistical error. For the s = 2 step scaling, the critical 
exponent of the beta function can be derived 7* = 0.31]'^Q'^g(stat.). This is also consistent 

with our main results with s = 1.5. 

Our result is consistent with -y*2-i°op ^ Q_3g ^j^^j loop(MS) ^ q 2g g^g estimated using 

2-loop and 4-loop (MS scheme) perturbation theory [51, 67, 68]. The result in the SF scheme 
is 7g^'^ = 0.13 ± 0.03 [8]. Our result provided a value larger than that in the SF scheme. We 
can conclude that both results are almost consistent with each other since the discrepancy 
of 7g is slight larger than l-u. 

Another scheme independent quantity which is interesting to observe is the mass anomalous 
dimension near the IRFP. That is the critical exponent for the relevant operator around the 
IRFP. We will report it in a forthcoming paper [69]. 

8. Summary 

We gave the explicit definition of the TPL renormalized coupling for SU(3) gauge theory and 
showed the running coupling constant in the case of quenched QCD and Nf = 12 theories. 
The definition is the extension of the SU(2) gauge theory, and we provided the perturbative 
calculation to define the coefficient in the case of SU(3). 

As a feasibility test for the SU(3) gauge theory, firstly we show the TPL running coupling 
constant in the case of quenched QCD. In the theory, there is confinement/deconfinement 
phase transition because of the finite volume effect and we study the behavior of the TPL 
renormalized coupling constant in the both phases. The TPL coupling have the remarkable 
property, that in the extremely low energy limit the coupling constant approaches to the 
constant (ffxPL ~ ™ case of SU(3) gauge theory), when the theory is in the confinement 
phase. From this analysis, the TPL coupling is found to be useful only in the deconfinement 
phase, so that we should take care of the phase in the parameter spaces. The running coupling 
constant is consistent with the perturbative result in high energy region, and it runs more 
slowly that that in the 1-loop perturbation in the low energy region. The running behavior 
is scheme dependent and is different from SF and Wilson loop schemes. 
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In the case of the theory coupled with fundamental fermions, we have to pay the atten- 
tion of the vacuum structure. We discussed the vacuum structure and the center symmetry 
breaking in our simulation setup using the semi-classical analysis, and generated the con- 
figurations in the true vacua. The vacuum structure depends on the boundary condition of 
the fermions, in the case of our definition, the configurations whose Polyakov loop in the 
untwisted direction has the nontrivial phase shows the minimum of the potential. 

Furthermore, we study the phase structure for Nf = 12 massive and massless fermion 
theories. There is a bulk phase transition near the massless region in /3 < 4.0. In such phase, 
the TPL coupling is not available since the theory shows the confinement behavior. 

Finally, we have found a solid evidence for the existence of an IRFP in SU(3) gauge 
theory with 12 massless Dirac fermions in the fundamental representation by using the TPL 
scheme. The coupling constant at the IRFP is 5xpl ~ 2.69. The stability of the fixed point 
is discussed, and we can conclude there is the IRFP after the systematic uncertainties are 
included. 
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A. Raw data of TPL coupling constant 
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Table Al L/a = 6, 8, and 10. (Data set A) 
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Table A2 L/a = 12, 16, and 20. (Data set A) (The data for /3 = 5.5, L/a = 12 becomes 
poor statistics since the bugged data was found after the simulations.) 
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Additional data set of TPL coupling constant for the local fit analysis 
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2.307(18) 


584700 


2.434(35) 


249700 


5.95 






2.317(64) 


74800 






5.90 










2.528(77) 


47300 


5.81 


2.248(11) 


530200 


2.471(22) 


415800 


2.605(37) 


216800 


5.80 






2.516(47) 


78200 






5.53 


2.408(11) 


718600 


2.676(29) 


471800 


2.784(46) 


177000 


5.36 


2.489(10) 


696700 


2.698(66) 


42900 


2.820(65) 


89400 



Table Bl L/a = 6, 8, and 10.(Data set B) 



31/41 





L/a -- 


= 12 


L/a -- 


= 16 




Q'-p p T 


# of Trj. 


Qnr p T 

i L^ L/ 


# of Trj. 


20.13 


0.3982(63) 


83000 


0.4138(75) 


79800 


17.55 


0.4662(68) 


85100 


0.4780(92) 


83600 


15.23 


0.588(12) 


75200 


0.566(14) 


80700 


13.85 


0.674(14) 


79100 


0.699(19) 


70800 


11.15 


0.914(13) 


173700 


0.962(30) 


72500 


9.42 


1.263(19) 


256200 


1.228(29) 


147600 


8.45 


1.470(18) 


397200 


1.520(25) 


368200 


7.82 


1.670(26) 


262300 


1.695(52) 


136300 


7.11 


1.966(30) 


256800 


1.996(40) 


244700 


6.76 


2.095(34) 


265700 


2.163(65) 


136400 


6.47 


2.278(29) 


330900 


2.391(50) 


273100 


6.12 


2.562(36) 


283200 


2.489(62) 


183700 


5.81 


2.723(45) 


269100 


2.729(79) 


186000 


5.53 


2.950(56) 


233600 


2.953(83) 


191200 


5.36 


3.030(50) 


272600 


3.06(10) 


187200 



Table B2 L/a = 12 and 16. (Data set B) 



C. Data set dependence, fit range dependence and the step scaling size 
dependence 

We show the two kinds of result for the growth rate from the global fit analysis and local 
fit analysis in Sec. 7.1 and 7.2 respectively. They are consistent with each other within 1-cr 
although they have the difference data sets, the different fit range and fitting function each 
other. Here we would like to show the each systematic uncertainties. 

The left top panel in Fig. CI shows the data set dependence. The central analysis includes 
both Data sets in Appendix A and B, while the blue result is obtained by only Data sets 
in Appendix A, which is the same data sets with the global fit analysis. The results are 
consistent with each other within l-cr. The right top panel in the Fig. CI shows the fit range 
dependence. The blue result is obtained by the data in the narrow /3 range; /3 < 7.0 for all 
lattice size. We can find the result is completely consistent with each other, and this local 
fit analysis is quite stable under the change of the fit range. Actually, we focus on the local 
P region where all data can be fitted by the linear or quadratic function in (3. Such results 
can be expected when the data can be fitted well. 

The bottom panel in Fig. CI shows the result of s = 2 step scaling. Although the step 
scaling function depends on the step scaling parameter, the comparison for the position u* 
must be independent of the step scaling parameter. We can find the position is consistent 
with that for s = 1.5, so that we can confidently conclude the existence of the fixed point 
and the interpolation works well in the s = 1.5 step scaling. 
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2 2.5 3 

U 



Fig. CI The comparison for several local fit result. The top left panel shows the data set 
dependence. The top right panel shows the fit range dependence. The bottom panel shows 
the s = 2 step scaling, and then the comparison for the position of u* shows the step scaling 
parameter dependence. 



D. The effect of the nonuniform data on the global fit 

The fixed point in this paper shows u* = 2.69, however, the paper [59] shows u* ~ 2.0. We 
consider that this difference comes from the mismatching the analysis method and the data 
sets quality in the paper [59]. In this appendix, we would like to consider the problem. 

The data in the paper [59] are strongly concentrated around /3 = 5.0 and they are a part 
of Data set A and all data of set B. In the paper [59], the authors carried out the global 
fit analysis. As we mentioned in our guiding principle point 1, when we use the global fit 
in the broad /3 region by using a single fitting function, ideally the data do not favor a 
specific region. Our analysis for the global fit used only the data A in which each data has 
a similar accuracy and the interval of /3 is chosen to realize the almost constant growth rate 
of the renormalized coupling on the lattice. On the other hand, the data set B is strongly 
concentrated around f3 = 5.0 - 6.0 and a part of the data has quite high statistics in the 
high /3 region and for the small L/a. In this appendix, we would like to derive the global 
step scaling function by using the data in the both sets A and B and discuss the effect of 
the prejudiced data. 

In Fig. Dl, the results by using the data set A and the ones by using both A and B are 
shown in the top and bottom panels respectively. Each procedure of the step scaling is the 
same with the Sec. 7.1 and Sec. 7.2 for the global and local fit analysis respectively. Each 
red result in fig. Dl denotes the global fit analysis and the blue one denotes the local fit 
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2-loop 
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Fig. Dl The data set dependence and the fit range dependence. The top panel is obtained 
by only the data set A, and the bottom one is by the data sets A and B. Each red data 
denotes the global fit analysis and the blue one shows the local fit analysis. We can find the 
bottom panel shows the large fit range dependence. 



analysis. The latter global fit result looks very similar with the result of the paper [59]. It 
shows the worse matching with the perturbative result, although still it is consistent with 
the perturbative prediction within the statistical error. The values u* becomes smaller than 
the other ones with more than 2-cr discrepancy. 

Now, we should consider which result is reliable. The global fit with the broad regime is 
always dangerous since the interpolated value includes non-vanishing contributions coming 
from the far region. To remove such effect we also carried out the local fit only with the 
focusing regime. The comparison the global fit result and the local fit result with both data 
A and B shows the larger fit range dependence rather than our central analysis in Fig. 17. 
The guiding principle point 1 must be important for such global fit analysis. That is the 
reason why we did not use the whole data of the data set B, when we would like to show 
the global behavior of the TPL running coupling constant. 

E. Comments on the estimation method of the discretization effect 

The discretization effect of the step scaling function T,{u,a/L; s) has two origins. The first 
one is a simple discretization effect of the renormalized coupling in the larger lattice size 
(sL/a). The second one comes from the tuned value of /? to reproduce the input quantity u 
in the smaller lattice (L/a). When we take the continuum limit, we fixed the physical box 
size L and the lattice spacing a (= and then the leading term of the former {0{{a/ sL)"^)) 
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is smaller than the later one (0((a/L)^)). So, it must be safe to avoid the interpolation of 
the small lattice sets if we consider the discretization effect seriously. 

In our simulation, we have L/a = 6, 8, 10, 12, 16 and 20. Then in the s = 1.5 step scaling we 
can take 4 data points to estimate the O(a^) effect with avoiding the small L/a interpolation 
on the other hand 5 = 2 step scaling has only 3 data points. One of the advantages to use 
s = 1.5 step scaling is that there is the finest lattice data(L/a = 12). Furthermore the chi 
square fit with only one degree of freedom is strongly disturbed the statistical fluctuation 
then we take the 4 points linear extrapolation as the central analysis in this paper. 

In the paper [59], they carried out the interpolation for L/a = 7 by using L/a = 6 and 8 
However, the interpolation is dangerous since it is the coarsest two lattices interpolation 
Actually, the raw data of the TPL (Figs. 13) shows the largest difference between L/a = 6 
and 8 in the low /3 region and that must induce a large uncertainty of the interpolation. 
Furthermore, the estimation of the systematic uncertainty between 4 data points linear 
extrapolation for L/a = 6, 7, 8, 10 — )• 12, 14, 16, 20 and 3 data points linear extrapolation for 
L/a = 7, 8, 10 —7- 14, 16, 20 is nonsense, since the data of L/a = 7 is generated by L/a = 6 
lattice data and thus the later extrapolation does not remove the effects of the coarsest 
lattice. 

F. Eigenvalue of the Dirac operator 

In this appendix, we report results for the measurement of the eigenvalue of the Dirac 
operator. We confirm two things from the quantity. At first, we show the global behavior of 
the eigenvalue. We find that the perturbative data is consistent with the tree level analysis in 
Sec. F.l. Furthermore, the /3 dependence of the data is smooth in the whole region and the 
data are inconsistent with zero even in the lowest f3 in our simulation parameter. These data 
also show that the theory is in the deconfinement and chiral symmetric phase. Secondly, we 
discuss the taste breaking in our simulation in Sec. F.2. In the case of the high /3 and the 
large lattice extent, the raw data of the low lying eigenvalues shows the degeneracy of the 
taste. In the strong coupling region, we take the continuum limit of the eigenvalue by using 
the TPL coupling constant as a reference of the same physics, and then we can show the 
effect of the taste breaking becomes mild in the continuum limit even in the strong coupling 
region. 

This section is a preliminary report for the study on the eigenvalues of the Dirac operator. 
We will report the detailed studies in the independent paper in near future [69]. 

F.l. Perturbative analysis 

Let us consider the massless staggered-Dirac operator D{x,y), 

D{x,y) = ^rif,{x) U^{x)6^j^[,^y - U^ix - fi)5^_^^y , (Fl) 

where r]^{x) is the staggered phase. This eigenvalues of D are pure imaginary since it is the 
anti-hermitian D^{x,y) = —D{x,y): 

D(x,y)^J)=zA«^f, (F2) 

where V'a^ denotes a Dirac fermion and (/) denotes the level of the eigenvalues. We define 
the lowest one as / = 1. The degree of freedom of one staggered-Dirac operator is 16, and 
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there are additional 3 color and 3 smell indices for our Dirac fermion (see: eq.(2.7)). The 
number of flavor (Nj = 12) is realized by (4 taste's) x (3 smell's) degrees of freedom and 
there are 3 flavor (=smell) staggered Dirac fermions. The operator D^D is hermitian and 
positive definite, and it can be decomposed to the operators on the even and odd sites. In 
this work, we measure the eigenvalues only positive and in the even-to-even sites, and then 
the degeneracy of one staggered-Dirac operator (4 taste's x 4 spinor's) becomes half. The 
remain degrees of freedom (3 color's x 3 smell's) can be counted as the unphysical twisted 
momenta in the twisted directions as explained in the Sec. 4. 



P=50.0 




1 1 I \ I \ I \ I \ I I 

20 40 60 80 100 

level (I) 



Fig. Fl The eigenvalues for /3 = 50, L/a = 6. The statistical error is the same size of the 
symbol. 



At the tree level, we can calculate the eigenvalue of the Dirac operator on the lattice: 

A^=4j;sin2^, (F3) 

where denotes the momentum of the fermion field for each direction. The leading order 
of 0{a) for low lying eigenvalues is proportional to the sum of k^. In our simulation there is 
twisted boundary condition and the non-trivial vacuum phase, and then the momentum is 
given by the eq. (4.3) as k^. The lowest momentum is given in the following case of 

i'^f!'^ ) = ^) °^ ^) ^ = X and y, 

= for both fi = z and t. 

The degree of degeneracy for these momentum combinations is 4, and the sum of k^ is given 
by y/TOir/L. The second lowest mode is given in the case of 

(—1, 2) or (0, 0) fov fj. = X or y, 
(—1, 1) or (0, 1) for fi = y or X (in same order), 
for both fj, = z and t 



(nr,o 



n' 



ph 
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The degree of degeneracy of that is 8, and the sum of the momentum is \f\^'K jL. The third 
lowest mode is also counted, it is given by 

(—1, 2) or (0, 0) for both [i = x and = y, 
for = z or t, 

— 1 for fj, = t or z (in same order). 

The number of degrees of degeneracy of that is 8, and the sum of the momentum is ^/22tt/L. 

Let us compare the measured value of the simulation with this tree level analysis. The 
total degree of the degeneracy should be multiplied by 8, since we measure the eigenvalue 
of the staggered fermion on even-to-even site as we explained. 

Figure Fl shows the first 100 eigenvalues in the /3 = 50, L/a = 6. There is a clear gap 
between / = 32 and / = 33 as we expected, although there is a large taste breaking since the 
lattice size is small. The ratio of the eigenvalue between first and 33rd levels is 1.328(2), and 
it is almost consistent with the tree level prediction ^J^■ On the other hand, we cannot see 
the second gap, which we expect to lie between 96th and 97th levels. The numerical value 
of the ratio of them is 1.110(1), and the value is completely consistent with y^f|- Although 
the second gap is not clear, the eigenvalue reproduces the tree level prediction. 



n 



ph 



n 



ph 



F.2. Global behavior of the eigenvalues and the taste breaking 

We also measure the eigenvalue for all lattice parameters in the data sets in Appendix A. We 
show some examples in Fig. F2, in which the error is estimated by the bootstrap analysis. 
The eigenvalues of high (3 and the large lattice size show the 4-fold degeneracy, although 
there is large taste breaking in the small lattice size or low /3 region. The data at the lowest 
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beta=6.0 




beta=16.0 


beta=5.0 




* beta=10.0 


* beta=4.5 





Fig. F2 The examples of the eigenvalues for several j3 and L/a. The statistical error is 
the same size of the symbol. 
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(3, (3 = 4.5 for L/a = 6, 12 and /3 = 5.5 for L/a = 20, show the inconsistency with zero, and 
the (3 dependence of the data at fixed lattice extent is smooth in whole /3 region. If we 
assume that the Banks-Casher relation, and the chiral symmetry is also preserving even in 
the lowest (3 in our simulation parameter. 

To see the taste breaking in our step scaling analysis in the strong coupling regime, we 
would like to show the eigenvalues in the continuum limit by taking the TPL renormalized 
coupling as a reference of the same physics: 



L- A 



(0 

cont 



limL- A(')(/3,L/a] 



(F4) 



Stpl {-^0 / L)=const 

where the (L • Xcont) is dimensionless quantity and the scale L is defined by the value of the 
renormalized coupling constant. The figures F3 show the (3 dependence of the eigenvalues 



1.2 



O 0.8 
cdl 0.6 - 

■-a 

0.' 



0.2 



level 1 



I I I I — r 



• L/a=6 
■ L/a=8 

L/a=10 
A L/a=12 
< L/a=16 

L/a=20 



I ' I ' I ' I 



• ■ 

• ■■■ ♦ ♦ ♦ * 



I I I I I I I 



♦ 

▲ 

T 



1.2 



0.4 



0.2 





level 33 


r' 1 ' 1 ' 1 


1 1 1 1 1 1 1 1 1^ 




• * 


• 


• 

• 

m 


• 

• 

//:""::: 

•■V ^^^^^ 


. ■ ■ ■ 


1 < 1 < 1 < 1 


< 1 < 1 < 1 < 1 < 1 



10 12 14 16 18 20 4 



10 12 14 16 

P 



18 20 



Fig. F3 The /3 dependence of the eigenvalues at the fixed L/a for the level 1 and 33. The 
statistical error is the same size of the symbol. 



for the level 1 and level 33 for each lattice size. We fit the data at fixed level and lattice size 
in terms of /3 by the fitting function, 

Af-l 

L-X^'\f3,L/a) = ^Ci/(3\ (F5) 

i=0 

where N is the number of the fitting parameter, and in practical we choose the best fit value 
of N for each / and L/a. Typically, N = A - 7 are employed in this analysis. 

In this analysis, the leading discretization error comes from the eigenvalue itself, which is 
proportional to L ■ A^'^ oc const. + O(a^) if the theory lives in the deconfinement phase. The 
other leading contribution O(a^) comes from the renormalized coupling, and the contribution 
comes via the tuned value of /3. We take the continuum limit for 5 data points, L/a = 
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8, 10, 12, 16 and 20, and they can be fitted well by the hnear function of (a/L)^ in whole 
region. To estimate the systematic uncertainty of this continuum extrapolation we also show 
the quadratic extrapolation of {a/L)"^ for 6 data points included the coarsest lattice L/a = Q. 



► level 8 

• level 7 
level 6 

■* level 5 

* level 4 

♦ level 3 
■ level 2 

• level 1 



I - 



0.01 0.02 0.03 

(a/L)' 

Fig. F4 The continuum extrapolation for the low lying eigenvalues for u = 2.5. Each full 
symbol denotes the eigenvalue at the interpolated (3 for each lattice size; L/a = 6, 8, 10, 12, 16 
and 20 from the right to left. The shadow symbol at (a/L)^ = denotes the extrapolated 
value by the linear extrapolation function of (o/L)^ for the finer 5 lattice data. The empty 
symbol at (a/L)^ = —0.001 shows the extrapolated value for the quadratic function of [ajVf' 
by using all 6 data points. 




u=2.5 



As an example, we take u = 2.5, which corresponds to the region for {P^L/a) = (5.4,6) 
- (6.4,20). The taste breaking of the raw data in these region is strong. The continuum 
extrapolation for u = 2.5 is shown in Fig. F4. In this plot, we show the eigenvalues at 
continuum limit for the low lying eigenvalues; level 1 ~ 8. The shadow symbol at (a/L)^ = 
denotes the extrapolated value for the linear extrapolation function of (a/L)^ for the finer 
5 lattice data. The empty symbol at (a/L)^ = —0.001 shows the extrapolated value for the 
quadratic function of (a/L)^ by using all 6 data points. The difference between these two 
kinds of the extrapolation can be identified as the systematic uncertainty. Including the 
systematic error, we find that the breaking of the level 1-4 becomes mild at the continuum 
limit even in the strong coupling regime. 

References 

[1] E. Bilgici, A. Flachi, E. Itou, M. Kurachi, C. -J. D. Lin, H. Matsufuru, H. Ohki and T. Onogi et ai, 

PoS LAT 2009, 063 (2009) 
[2] E. Itou, T. Aoyama, M. Kurachi, C. -J. D. Lin, H. Matsufuru, H. Ohki, T. Onogi and E. Shintani et 

al, PoS LATTICE 2010, 054 (2010) 
[3] T. Aoyama, H. Ikeda, E. Itou, M. Kurachi, C. -J. D. Lin, H. Matsufuru, K. Ogawa and H. Ohki et al., 

arXiv: 1109.5806 [hep-lat]. 
[4] Y. Iwasaki, K. Kanaya, S. Kaya, S. Sakai and T. Yoshie, Phys. Rev. D 69, 014507 (2004) 



39/41 



S. Catterall and F. Sannino, Phys. Rev. D 76, 034504 (2007) 

S. Catterall, J. Giedt, F. Sannino and J. Schneible, JHEP 0811, 009 (2008) 

S. Catterall, J. Giedt, F. Sannino and J. Schneible, arXiv:0910.4387 [hep-lat]. 

T. Appelquist, G. T. Fleming and E. T. Neil, Phys. Rev. Lett. 100 (2008) 171607, [Erratum-ibid. 102 
(2009) 149902], Phys. Rev. D 79 (2009) 076010 

T. Appelquist et al. [LSD Collaboration], Phys. Rev. Lett. 104, 071601 (2010) 

T. Appelquist, G. T. Fleming, M. F. Lin, E. T. Neil and D. A. Schaich, Phys. Rev. D 84 (2011) 054501 
T. Appelquist, R. C. Brower, M. I. Buchoff, M. Cheng, S. D. Cohen, G. T. Fleming, J. Kiskis and 
M. Lin et al, arXiv: 1204.6000 [hep-ph]. 

Z. Fodor, K. Holland, J. Kuti, D. Nogradi and C. Schroeder, Phys. Lett. B681, 353-361 (2009). 
Z. Fodor, K. Holland, J. Kuti, D. Nogradi and C. Schroeder, PoS LAT 2009, 058 (2009) 
Z. Fodor, K. Holland, J. Kuti, D. Nogradi and C. Schroeder, JHEP 0911, 103 (2009) 
Z. Fodor et al, Phys. Lett. B 703 (2011) 348 

A. Deuzeman, M. P. Lombardo and E. Pallante, Phys. Lett. B 670, 41 (2008) 

A. Deuzeman, M. P. Lombardo and E. Pallante, Phys. Rev. D 82, 074503 (2010) 

K. Miura, M. P. Lombardo and E. Pallante, Phys. Lett. B 710, 676 (2012) 

A. Deuzeman, M. P. Lombardo, T. N. da Silva and E. Pallante, arXiv; 1209. 5720 [hep-lat]. 

P. de Forcrand, S. Kim and W. linger, arXiv:1208.2148 [hep-lat]. 

X. -Y. Jin and R. D. Mawhinney, PoS LATTICE 2011, 066 (2011) 

A. Hasenfratz, Phys. Rev. D 82, 014506 (2010) 

A. Hasenfratz, arXiv:1106.5293 [hep-lat]. 

A. Cheng, A. Hasenfratz and D. Schaich, Phys. Rev. D 85, 094509 (2012) 
T. DeCrand and A. Hasenfratz, Phys. Rev. D 80, 034506 (2009) 
T. DeCrand, Y. Shamh and B. Svetitsky, Phys. Rev. D 82, 054503 (2010) 
T. DeCrand, Phys. Rev. D 84, 116901 (2011) 

T. DeCrand, Y. Shamir and B. Svetitsky, Phys. Rev. D 85, 074506 (2012) 
T. DeCrand, Y. Shamir and B. Svetitsky, arXiv:1201.0935 [hep-lat]. 

Y. Aoki, T. Aoyama, M. Kurachi, T. Maskawa, K. -i. Nagai, H. Ohki, A. Shibata and K. Yamawaki et 
al, Phys. Rev. D 86, 054506 (2012) 

L. Del Debbio, A. Patella and C. Pica, Phys. Rev. D 81, 094503 (2010) 
L. Del Debbio, B. Lucini, A. Patella, C. Pica and A. Rago, Phys. Rev. D 80, 074507 (2009) 
L. Del Debbio, B. Lucini, A. Patella, C. Pica and A. Rago, Phys. Rev. D 82, 014510 (2010) 
F. Bursa, L. Del Debbio, L. Keegan, C. Pica and T. Pickup, Phys. Rev. D 81, 014505 (2010) 
L. Del Debbio and R. Zwicky, Phys. Rev. D 82, 014502 (2010) 

F. Bursa, L. Del Debbio, L. Keegan, C. Pica and T. Pickup, Phys. Lett. B 696, 374 (2011) 
L. Del Debbio and R. Zwicky, Phys. Lett. B 700, 217 (2011) 

F. Bursa, L. Del Debbio, D. Henty, E. Kerrane, B. Lucini, A. Patella, C. Pica and T. Pickup et al, 
Phys. Rev. D 84, 034506 (2011) 

A. Patella, Phys. Rev. D 86, 025006 (2012) 
M. I. Buchoff, Phys. Rev. D 85, 074503 (2012) 

M. Hayakawa, K. I. Ishikawa, Y. Osaki, S. Takeda, S. Uno and N. Yamada, Phys. Rev. D 83 (2011) 
074509 

T. DeCrand, Y. Shannr and B. Svetitsky, Phys. Rev. D 82, 054503 (2010) 
J. Ciedt and E. Weinberg, Phys. Rev. D 84, 074501 (2011) 
J. Ciedt and E. Weinberg, Phys. Rev. D 85, 097503 (2012) 

A. J. Hietanen, J. Rantaharju, K. Rummukainen and K. Tuominen, JHEP 0905, 025 (2009) 

A. J. Hietanen, K. Rummukainen and K. Tuominen, Phys. Rev. D 80, 094504 (2009) 

T. Karavirta, A. Mykkanen, J. Rantaharju, K. Rummukainen and K. Tuominen, JHEP 1106, 061 

(2011) 

T. Karavirta, J. Rantaharju, K. Rummukainen and K. Tuominen, JHEP 1205, 003 (2012) 
T. Karavirta, K. Tuominen and K. Rummukainen, Phys. Rev. D 85, 054506 (2012) 
W. E. CasweU, Phys. Rev. Lett. 33, 244 (1974). 

J. A. M. Vermaseren, S. A. Larin and T. van Ritbergen, Phys. Lett. B 405 (1997) 327 

T. Banks and A. Zaks, Nucl. Phys. B 196, 189 (1982). 

M. Luscher, P. Weisz and U. Wolff, Nucl. Phys. B 359 (1991) 221. 

M. Luscher, R. Narayanan, P. Weisz and U. Wolff, Nucl. Phys. B 384 (1992) 168 

M. Luscher, R. Narayanan, R. Sommer, U. Wolff and P. Weisz, Nucl. Phys. Proc. Suppl. 30, 139 (1993). 

G. M. de Divitiis, R. Frezzotti, M. Guagnelh and R. Petronzio, Nucl. Phys. B 422 (1994) 382 
G. 't Hooft, Nucl. Phys. B 153, 141 (1979). 

M. Liischer and P. Weisz, Nucl. Phys. B 266 (1986) 309 

C. -J. D. Lin, K. Ogawa, H. Ohki and E. Shintani, JHEP 1208, 096 (2012) 



40/41 



[60] G. M. de Divitiis, R. Frezzotti, M. Guagnelli and R. Petronzio, Nucl. Phys. B 433 (1995) 390 Nucl. 

Phys. B 437 (1995) 447 
[61] G. Parisi, Published in Cargese Summer Inst. 1983:0531 

[62] H. D. Trottier, N. H. Shakespeare, G. P. Lepage and P. B. Mackenzie, Phys. Rev. D 65 (2002) 094502 
[63] G. de Divitiis et al. [Alpha Collaboration], Nucl. Phys. B 437, 447 (1995) 

[64] S. Capitani, M. Luscher, R. Sommer and H. Wittig [ALPHA Collaboration], Nucl. Phys. B 544 (1999) 
669 

[65] E. Bilgici, A. Flachi, E. Itou, M. Kurachi, C. -J D. Lin, H. Matsufuru, H. Ohki and T. Onogi et al, 

Phys. Rev. D 80, 034507 (2009) 
[66] E. T. Tomboulis, arXiv:1211.4842 [hep-lat]. 
[67] T. A. Ryttov and R. Shrock, Phys. Rev. D 83 (2011) 056011 

[68] M. Mojaza, C. Pica and F. Sannino, Phys. Rev. D 82, 116009 (2010), C. Pica and F. Sannino, Phys. 

Rev. D 83, 035013 (2011) 
[69] E. Itou and A. Tomiya, in preparation 



41/41 



